Phase diagram of silicon from atomistic simulations 
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In this letter we present a calculation of the temperature-pressure phase diagram of Si in a 
range of pressures covering from —5 to 20 GPa and temperatures up to the melting point. The 
phase boundaries and triple points between the diamond, liquid, /3-Sn and Si34 clathrate phases 
are reported. We have employed efficient simulation techniques to calculate free energies and to 
numerically integrate the Clausius-Clapeyron equation, combined with a tight binding model capable 
of an accuracy comparable to that of first-principles methods. The resulting phase diagram agrees 
well with the available experimental data. 

PACS numbers: 64.30.+t, 64.70.-p, 65.40.-b 
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Temperature-pressure phase diagrams charter the re- 
gions of stability of the different allotropes of a material. 
The confection of phase diagrams has been a long stand- 
ing objective of experimental physics, chemistry and ma- 
terials science. However, to date, the phase diagrams 
of most materials remain relatively unknown beyond the 
domain of normal conditions, because of the technical 
challenge of performing accurate phase behavior stud- 
ies in conditions of extreme temperatures and/or pres- 
sures. Reliable first principles electronic structure cal- 
culations have the potential to be of great assistance 
in this problem, and indeed they have proved their value 
with impressive demonstrations of their capabilities, such 
as the calculation of the melting curve of iron down to 
the pressure regime of the Earth's core 0], that of alu- 
minum or or that of hydrogen in a similar range 
of pressures Q. But such calculations, which employ 
either free energy evaluation techniques like thermody- 
namic integration, or directly address phase coexistence 
by explicitly simulating the interface, are computation- 
ally demanding, and by no means routine. The two- 
phase method, in particular, requires large simulation 
cells where the two phases can be monitored in coex- 
istence, and it is only directly applicable to solid-liquid 
equilibria. Nevertheless, in recent years several simula- 
tion techniques have been developed which now make free 
energy calculations 0.01 an d phase boundary determina- 
tion |2j much more accessible. In this letter we demon- 
strate the potential of these novel techniques by using 
them to obtain, entirely from atomistic simulations, the 
phase diagram of Si in a wide range of temperatures and 
pressures. 

In spite of being one of the most extensively stud- 
ied materials, the phase diagram of Si is not accurately 
known. As many as eleven phases other than diamond 
have been identified at high pressures ||, with at least 
six of them being thermodynamically stable in some 
temperature-pressure domain. It is now also clearly es- 
tablished that other phases become stable at negative 
pressures: the so called clathrate phases .10], of which 



Si34 (also known as Sii36) has interesting semiconduct- 
ing properties ^lj of its own. In Si, as in many other 
materials, computer simulations have been a great aid 
in identifying, and sometimes even predicting |l2l 
the occurrence of certain phases, but they have not been 
much used to help establishing the limits of stability of 
the different phases except at zero temperature or 
zero pressure In this letter we calculate the 

phase diagram of Si entirely from atomistic simulations, 
for temperatures ranging between and 1700 K, and for 
pressures in the range —5 to 20 GPa. We have considered 
four different phases in this temperature-pressure region, 
namely the diamond structure (Si I), the /3-Sn phase (II), 
S134 (C) and the liquid phase (L), and provide the five 
corresponding coexistence curves between these, as well 
as estimates for the location of the two triple points to 
be found in this area of the phase diagram. 

In a recent study 01 > we have shown that certain tight 
binding 0] models, such as those of Kwon et al. and 
Lenosky et al. [2jj, are capable of providing very accu- 
rate descriptions not only of the structural properties of 
Si, but more importantly of the thermal properties too. 
In fact, these two models predict a melting temperature 
at zero pressure which is in better agreement with the 
experimental value than that provided by first principles 
calculations 0, Q|. The deciding factor between the 
Kwon and Lenosky models is that only the latter cor- 
rectly predicts a pressure-induced transition from the I 
to II phase [2(|- All calculations reported below were 
carried out using supercells containing 128 Si atoms, ex- 
cept in the case of Si C, where a supercell of 136 atoms 
was imposed by the structure. Four special k-points [2l| 
were used to sample the Brillouin zone and provided 
a sufficient degree of convergence even for the metallic 
phases (II and L). All calculations were performed with 
the Trocadero code |2^ . 

Let us now briefly describe the procedure adopted for 
determining the phase diagram. Firstly, for all pairs of 
phases for which a phase boundary is sought, a coexis- 
tence point along the boundary must be found, i.e. a tern- 
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Phase boundary 


Pressure 


Temperature 


dT/dP 




(GPa) 


(K) 


(K GPa" 1 ) 


I-L 





1551 ± 66 


-16 






(1687°) 


(-33°) 


II-L 


10 


1230 ± 25 


40 


C-L 





1424 ± 57 


-73 






(1473 6 ) 




I-II 


15.5 ±0.4 


500 


-200 




(10.4 - 12.4 C ) 


(573 c ) 




I-C 


-2.47 ±0.03 


500 


« 10 4 



a Reference^, b reference^, c reference !251 

TABLE I: Coexistence points obtained between the different 
phases considered in this study. Experimental data is given in 
parenthesis when available. These points were later used as 
starting initial conditions for dynamical Clausius-Clapeyron 
integration. The local slope of the corresponding coexistence 
line is also given. 



perature and pressure for which the Gibbs free energies 
of the two phases in question are equal. If the coexistence 
line has a small (in absolute value) pressure derivative, 
in order to locate the coexistence point it is better to fix 
the pressure of both phases at some convenient value, P, 
and to calculate the Gibbs free energy of each phase in a 
temperature interval bracketing the coexistence temper- 
ature, T c , at fixed pressure. This can be done employing 
the reversible scaling technique of de Koning et al. [6|. 
If, on the other hand, the phase boundary is expected to 
have a large (in absolute value) pressure derivative, as is 
common in solid-solid coexistence lines, it is more con- 
venient to fix the temperature and to monitor the Gibbs 
free energy of each phase as a function of pressure, which 
can be done with the adiabatic switching technique 
Once a coexistence point has been thus located, the rest 
of the phase boundary is obtained by solving numerically 
the Clausius-Clapeyron equation: 



dT c _ AV 

~dp ~ Tc ~Kh' 



(1) 



where T c is the coexistence temperature at pressure P, 
and AV and AH are the difference of volumes and en- 
thalpies of the two phases, respectively. This was done 
using the dynamical scheme of de Koning and cowork- 
ers 7]. All the above techniques require to simulate 
the system under isothermal-isobaric conditions, and this 
was done as described by Hernandez j^. Following 
the scheme outlined above, we proceeded to locate ini- 
tial coexistence points along the phase boundaries be- 
tween the diamond-liquid (I— L), /3-Sn-liquid (II— L), 
clathrate-liquid (C— L), diamond-clathrate (I— C) and 
diamond-/3-Sn (I— II) phases. Further calculations were 
carried out at the found coexistence points in order to 
quantify the errors incurred in our estimations of those 



conditions. For coexistence temperatures at fixed pres- 
sure the error can be estimated from ST C w SG/AS [24|. 
where SG is the error in the Gibbs free energy difference, 
and AS is the entropy change. For coexistence pressures 
at fixed temperature, the error can be estimated from 
SP 6G/AV, where Ay is the change of volume. As 
well as error estimates, these calculations allowed us to 
obtain the values of the pressure derivatives of each phase 
boundary at the coexistence point from Eq. 0). Table JIJ 
lists the locations of the different coexistence points, to- 
gether with their estimated errors and the local pressure 
derivatives of the corresponding coexistence lines. 

Our calculated zero-pressure melting point for Si-I 
(1551 K) is in reasonably good agreement with the ac- 
cepted experimental value of 1687 K @|. Although the 
difference between these values may appear to be large, 
we note that density functional theory (DFT) calcula- 
tions using the local density approximation (LDA) for 
the exchange-correlation energy predict values in the 
range 1300-1350 K [H El, and 1492 ± 50 K when a 
generalized-gradient approximation (GGA) is used in- 
stead [lij . Comparison with our own previous estima- 
tion using the same model ^3 gives a difference of 30 K, 
within our error bars for the melting temperature, al- 
though we regard the present value as more accurate. In 
agreement with experiments and with DFT calculations 
we obtain a negative slope for the phase boundary. The 
melting of (metastable) Si-II at 10 GPa is found to occur 
at 1230 K. This temperature is already higher than the 
experimental estimate of the I— II— L triple point temper- 
ature 25], but is consistent with our own prediction for 
this triple point (see below). Like the I phase, C has a 
melting line with a negative slope. For this phase we pre- 
dict a melting temperature at zero pressure of 1424 K, 
which is only 50 K below the experimentally measured 
value at this pressure 0, and also in good agreement 
with the calculated value of Wilson & McMillan As 
for the coexistence pressure between Si-I and Si-II at 
500 K, 15.5 GPa, it is larger than the experimental val- 
ues, 10.4-12.4 GPa at 573 K [2j|, but we note again that 
DFT-LDA calculations predict a value of 8 GPa Q at 
zero temperature, while recent quantum Monte Carlo cal- 
culations place it at 16.5 GPa |2jj. At the same temper- 
ature, the I— C coexistence point is found at —2.47 GPa, 
in good agreement with experimental measurements 
and with both empirical potential simulations |27l ] and 
first principles calculations |2*fj| . 

Taking as starting conditions the coexistence points 
thus located, we then proceeded to run dynamical 
Clausius-Clapeyron integration calculations, thereby 
obtaining the sought phase boundaries. Fig. JIJ shows 
our calculated phase diagram for Si, and constitutes the 
central result of this work. For comparison, Fig. also 
shows a summary of recent experimental observations 
from Voronin et al. .25], Hu et al. jH and McMillan 0. 
As can be seen, the calculated phase diagram captures 
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all the main features of the experimental data with sur- 
prising fidelity. Nevertheless, there are differences in the 
details which are mostly attributable to minor shortcom- 
ings in the Lenosky tight-binding model. 

The I— L coexistence line has, as previously mentioned, 
a negative pressure derivative, which increases slightly 
toward larger negative values as the pressure is raised. 
This is in agreement with experimental observations by 
Voronin et al. [2^ that suggest this behavior, although 
the pressure derivative predicted by the Lenosky model 
for this phase boundary is too small. An independent 
error estimate of the melting line away from the start- 
ing point of the dynamical Clausius-Clapeyron integra- 
tion calculation was obtained at a point of coordinates 
T = f385 K and P = 8 GPa. We found an uncer- 
tainty of 95 K, which is not significantly worse than 
that of the zero- pressure melting point. The I— L melt- 
ing line meets the I— II phase boundary at a temperature 
T w f290 K and P « 10.9 GPa, according to our re- 
sults. The experimental coordinates of this triple point 
are not accurately known, though a recent estimate by 
Voronin et al. 25] puts it at T = 1003 ± 20 K and 
P = 10.5 ± 0.2 GPa. Compared with this best exper- 
imental estimate, our triple point temperature is some- 
what too high (by nearly 300 K), consistent with the fact 
that dT c /dP for the I— L melting line is too small in ab- 
solute value. However, we emphasize that the estimate of 
Voronin et al. is a lower bound; if one assumes that the I— 
L value of dT/dP remains constant and equal to its zero 
pressure value, then at 10.5 GPa the triple point tem- 
perature should be 1340 K, which is much closer to our 
figure. Thus, it is very likely that our triple point temper- 
ature and that of Voronin et al. provide upper and lower 
bounds respectively for the true value. The agreement in 
the value of the pressure coordinate is much better (in 
fact, within our error estimate for the I— II phase bound- 
ary), but must be understood as somewhat fortuitous, 
resulting from error cancellation between a I— II coexis- 
tence pressure at 500 K which is slightly too high, and 
a value of dT c /dP for the I— II phase boundary which is 
most likely smaller than the experimental one. Never- 
theless, it should be pointed out that, to our knowledge, 
this is the first prediction from atomistic simulations of 
the location of the of the I— II— L triple point of Si. We 
have also calculated the II L phase boundary, starting 
from a temperature and pressure where both the II and 
L phases are metastable. Indeed, the II L coexistence 
line crosses the I— II and I— L boundaries very close to the 
point where the I— L and I— II boundaries cross, a good 
indication of the internal consistency of our calculations. 

The I— II phase boundary has, as expected, a large 
pressure derivative, of negative sign. Experimental ob- 
servations [25| seem to agree with this finding, although 
data is only available at large temperature intervals, and 
it is not possible at present to compare with an exper- 
imental value of the pressure derivative. In any case, 
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FIG. 1: Silicon phase diagram. The continuous and dashed 
lines indicate the calculated phase diagram. Dashed curves in- 
dicate phase boundaries in regions where the separated phases 
are metastable, while continuous black curves separate ther- 
modynamically stable phases; uncertainty bounds estimated 
at specific points of the phase diagram (marked by filled cir- 
cles) are provided by the error bars. For comparison pur- 
poses, a schematic phase diagram summarizing the exper- 
imental data is shown with dotted lines, and experimental 
data at specific temperatures and pressures is shown in the 
form of empty symbols. The asterisk corresponds to the zero- 
pressure melting point of phase I, 1687 K 0; the circle is 
the zero-pressure melting point of the (metastable) C phase, 
at 1473 K 0; the diamond is the I— II— L triple point, with 
estimated coordinates of 1003 ±20 K and 10.5 ±0.2 GPa |2j|; 
empty squares and triangles indicate the pressures at which 
the II phase was first observed and where the I phase ceased 
to be detected, respectively, in the experiments of Voronin 
et al. [25| : left and right pointing triangles give the same in- 
formation as obtained by Hu et al. |30||: finally, the down- 
ward pointing triangle is the estimated I— C— L triple point, 
at 1710 K and -2.5 GPa H. 



it can be concluded that, as it happened with the I— L 
melting curve, the pressure derivative of the I— II coex- 
istence line is probably not as large as the experimental 
one. This can be seen from the fact that both coexis- 
tence pressures at T = 500 K and in the limit T-*0K 
are several GPa higher than the experimental coexistence 
pressures at similar temperatures, while our pressure co- 
ordinate for the I— II— L triple point agrees quite well with 
the experimental value. It is worth noting that the coex- 
istence pressure at T — ► K deduced from the calculated 
I— II phase boundary (18.1 GPa) is very close to the value 
obtained by calculating the enthalpies of both phases at 
T = K (18.5 GPa), which serves as another internal 
consistency check for our calculations. To further vali- 
date our results, we performed an error estimation of the 
coexistence pressure at 1000 K. Our results indicate that 
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the coexistence pressure predicted by the phase bound- 
ary at that temperature (12.6 GPa) is accurate to within 
0.4 GPa. 

The I— C phase boundary we have obtained is nearly 
a vertical line at P « —2.5 GPa. At 500 K our calcu- 
lated coexistence point occurred at a pressure of — 2.47± 
0.02 GPa [see Table |JTJ]. At 1000 K, we calculated a sec- 
ond coexistence point to double check the phase bound- 
ary calculation, obtaining a pressure of —2.4 ± 0.4 GPa, 
to be compared with a value of —2.5 GPa according to 
our calculated phase boundary. The pressure derivative 
of this phase boundary is so large (« 10 4 K/GPa in ab- 
solute value) as to make it virtually impossible to pre- 
dict its sign with any accuracy. This is a manifestation 
of the fact that at coexistence, not only the Gibbs free 
energies of the two phases are the same, but also their 
enthalpies are nearly equal. This in turn implies that the 
entropies of the I and C phases are very similar at coex- 
istence conditions. The I— C and I— L boundaries cross at 
T « 1576 K and P w -2.48 GPa, marking within our er- 
ror bars the coordinates of the I— C— L triple point. Wil- 
son & McMillan j2]j have estimated the location of this 
triple point to be 1750 K and —1.5 GPa, from two-phase 
coexistence calculations employing the Stillinger- Weber 
potential j^. Experimental estimates @ suggest that 
the triple point may actually be closer to 1710 K and 
—2.5 GPa. Again, we note that our predicted tempera- 
ture coordinate is too low by about 100-150 K, consistent 
with our underestimation of the I melting temperature 
at zero pressure, while the pressure coordinate is closer 
to the experimental estimation. The C— L melting line is 
almost straight, with a negative slope of —73 K/GPa; it 
crosses the I— L melting line slightly to the right of the 
point where the latter is crossed by the I— C boundary, 
but their separation is within our estimated error bars. 

Thus it is seen that, in spite of its simplicity and semi- 
empirical nature, the Lenosky model provides a fairly 
good description of the phase diagram of Si, being prob- 
ably as accurate as could be expected of first principles 
calculations. We have obtained five phase boundaries 
and two triple points between four phases of the silicon 
phase diagram, in reasonable agreement with the known 
experimental data. The simulation techniques employed 
in this study to calculate free energies and to obtain coex- 
istence curves are straight forward and efficient, and can 
equally well be used in combination with first principles 
methods. It can be concluded, then, that the combina- 
tion of techniques used here brings about the possibility 
of obtaining entire phase diagrams of complex materials 
completely ab initio. 
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